The tail of the contact force distribution in static granular materials 
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We numerically study the distribution P{f) of contact forces in frictionless bead packs, by aver- 
aging over the ensemble of all possible force network configurations. We resort to umbrella sampling 
to resolve the asymptotic decay of P[f) for large /, and determine P{f) down to values of order 
10~*^ for ordered and disordered systems in two and three dimensions. Our findings unambiguously 
show that, in the ensemble approach, the force distributions decay much faster than exponentially: 
P{f ) ~ exp(— with a~2.0 for 2D systems, and a~1.7 for 3D systems. 

PACS numbers: 45.70.Cc, 05.40.-a, 46.65.+g 



The contact forces inside a static packing of grains are 
organized into highly heterogeneous force networks, and 
can be characterized by the probability density of contact 
forces P(/) Q. Such force statistics were first studied 
in a series of experiments that measured forces through 
imprints on carbon paper at the boundaries of a granu- 
lar assembly. Unexpectedly, the obtained P{f) displayed 
an exponential rather than a Gaussian decay for large 
forces 0. After these initial findings, other experimental 
techniques have revealed similarly exponentially decay- 
ing distributions of the boundary forces Q] • 

The first model that captured this exponential decay 
was the pioneering g-model, where scalar forces are bal- 
anced on a regular grid Later studies found, how- 
ever, that the nature of the tail of P{J) depends on the 
details of the stochastic rules for the force transmission 
in this model and need not be exponential Other 
explanations for the exponential tail hinge on "entropy 
maximization" , or closely related, on an analogy with 
the Boltzmann distribution [^, . The essence of the lat- 
ter argument is that a uniform sampling of forces that (i) 
are all positive (corresponding to the repulsive nature of 
contact forces), and (ii) add up to a constant value (set 
by the requirement that the overall pressure is constant), 
strongly resembles the microcanonical ensemble, in which 
configurations are flatly sampled under the constraint of 
fixed total energy. 

As it is difficult to experimentally access contact forces 
inside the packing, many direct numerical simulations of 
P{f) have been undertaken [13, El- While numerous of 
these studies claim to find an exponential tail as well, 
the evidence is less convincing than for the carbon paper 
experiments: apart from [lo| . nearly all numerical force 
probabilities bend down on a logarithmic plot, suggest- 
ing a faster than exponential decay [llj . In addition, new 
experimental techniques using photoelastic particles 12 1 
or emulsions [3, [3], have produced bulk measurements, 
and these also reveal a much faster than exponential de- 



cay for P{f), consistent with a Gaussian tail. 

These contradicting findings completely reopen the 
discussion on the tail of P{f). The presently available 
data for P{f) have been obtained from a wide variety 
of systems and models, and parameters such as dimen- 
sionality, hardness of grains, bulk vs boundary measure- 
ments, may ultimately all play a role in determining the 
asymptotics of P{f). In addition, in many cases, the true 
asymptotic nature of P{f) is hard to probe, since obtain- 
ing reliable data for forces much larger than the average 
value (/) remains a challenge. 

In this paper we will probe the tail of P{f) in the force 



network ensemble [H |T3, 0, El . We numerically 
resolve the probability for large forces using the technique 
of umbrella sampling [2l[ , which yields accurate statistics 
for P{f) for relative probabilities down to 10"''^ and / up 
to / = 15(/). This high accuracy is crucial for excluding 
any cross-over effects and allows to unambiguously 
identify the behaviour for f ^ (/) . We study frictionless 
systems in two and three dimensions, both with ordered 
and disordered contact networks, and also explore the 
effect of system size and contact number. 

For all these systems, we have found that the ensemble 
yields a much faster than exponentially decaying force 
distribution. The dimensionality of the system is crucial, 
while other factors hardly affect the asymptotics: P{f) 
decays as exp(— /"), with Q! = 2.0±0.1 in two dimensions. 
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while in three dimensions a = 1.7 ± 0.1 

Force network ensemble and umbrella sampling — 
The ensemble approach to force networks is inspired by 
the proposal of Edwards to assign an equal probability 
to all "blocked" states, i.e. states that are at mechani- 
cal equilibrium j23l |. By limiting the Edwards ensemble 
to a sin gle contact geometry of a packing of frictionless 
spheres 2^, where the contact forces are the remaining 
degrees of freedom and all allowed force-configurations 
are sampled with equal weight, we obtain the force net- 
work ensemble. We restrict ourselves to spherical par- 
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tides with frictionless, repulsive contacts, so that every 
contact force fi corresponds to one scalar degree of free- 
dom. Furthermore, we require all fi > due to the 
repulsive nature of the contacts. As the equations of 
mechanical equilibrium are linear in the contact forces, 
one can cast the solutions / = (/i, /2, • • •) in the form 
f = fo + '^if CkVk- The solution space is spanned by 
the vectors Vk and /o, and can be sampled through the 
coefficients Ck - for details we refer to Refs. [15, 16, 1^. 
For a hexagonal packing (two dimensional), these vectors 
are easily constructed using so-called "wheel moves" [31 , 
but for other packings we have obtained Vk and /o from 
a simulated annealing procedure 15]. Ensemble averages 
using a uniform measure in this force space then become 



(1) 



(q) =n-' dcq , n= dc, 



where the integral runs over the coefficients Ck limited to 
the convex subspace C for which all fi>0 [11]. 

To obtain accurate statistics for large forces we per- 
form umbrella sampling. The idea is to bias the nu- 
merical sampling towards solutions with large forces, us- 
ing a Monte Carlo technique with a modified measure 
dcp{c)/VL, and then correct for this bias when perform- 
ing the averages, {q) = {q / p) umhrdu, since 



(q) = n- 



dcp{c) . , 



(2) 



Defining /max as the largest force for a given c, we have 
used a measure p{c) oc e'^^'''"""'^ where W is chosen such 
that the probability of /max in the modified ensemble 
is approximately flat in the range (/) < /max < 15(/}. 
We have verified that this procedure exactly reproduces 
P(/) in the range accessible by the conventional unbiased 
sampling [l^. However, forces of the order of 15(/) are 
now sampled only lO'* times less frequently than around 
(/), even though their relative probability is about lO""*^, 
leading to the spectacular improvement in the numerical 
accuracy [25| . 

Hexagonal packings in two dimensions — A well stud- 
ied geometry for which the force network ensemble yields 
nontrivial results is when all particles are of equal size 



and form a hexagonal lattice [15|, H, 1^. The umbrella 
sampling allows us, for the first time, to access the statis- 
tics beyond / = 5(/}. Figure [Tfa) shows that P{f) de- 
cays much faster than exponentially, and that effects of 
the finite size of the system are weak. Figure [TJb) illus- 
trates that for increasingly large systems, P{f) rapidly 
converges to an asymptotic form which is characterized 
by a purely Gaussian decay. This can also be seen in the 
inset of Fig. 1(b), where we exploit the fact that we have 
access to P{f) over more than forty decades: Assuming 
that for large /, P{f) ~ exp(— (//A)"), one can infer the 
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FIG. 1: Force probabilities in two dimensional, hexagonal 
packings with periodic boundary conditions and increasing 
system size as indicated by the number of particles A^. (a) 
P(/) decays much faster than exponentially, and rapidly con- 
verges to its asymptotic form with A'^. The inset illustrates 
that system size effects are hardly visible for P(/) down to 
10~®. (b) logP vs becomes a perfectly straight line for 
large systems, indicating that the tail of P(/) is well described 
by a Gaussian decay ~ exp(— /^) (dashed line). The inset 
shows that on a triple- log plot, the asymptotic decay attains 
a slope close to 2, confirming the Gaussian tail (see text). 
Curves are offset for clarity, and lines are guides to the eye. 



exponent a from the asymptotic slope of a triple-log plot 
in which log(— log P) is plotted as function of log / (base 
10) In Fig. [He) we find a to be 2.0 ± 0.1, confirm- 
ing that the tail of P{f) is well described by a Gaussian 
decay 

Disordered packings in two dimensions — To inves- 
tigate the effect of packing disorder and coordination 
number z, we have created packings from molecular dy- 
namics simulations of soft particles in periodic boundary 
conditions (see l^, 17|). The coordination number z is 
controlled by the pressure in the simulations. Once a 
packing is obtained, its geometry is kept fixed, and we 



subsequently explore the ensemble of force networks for 
these packings. 

For all 2D disordered packings, P{f) decays much 
faster than exponentially, as shown in Fig. [21 Comparing 
the ordered hexagonal packings to a disordered system 
with equal coordination number, z = 6, we find nearly 
indistinguishable P[f) (inset Fig. 2(a)). This suggests 
that the packing (dis)order is not important for P{f). 

The contact number influences the asymptotics of 
P{f )'- a lower z leads to a faster decay (Fig. [Ha-b)), 
although in the restricted range / < 5(/), the force dis- 
tribution appears very close to Gaussian for all z (inset 
Fig. [H^b)). For the lowest z in particular, this tendency 
is cut off at large /, which can be clearly seen in the 
triple-log plot Fig. [Ifc), where all curves tend towards 
a well-defined slope a ~ 2 for intermediate / but cross 
over to a much faster decay for large /. We suggest that 
this is a finite size effect, which is most severe when z 
approaches the isostatic point (z — 4), where there are 
less and less degrees of freedom available [illilj. Indeed, 
data for z = 4.5 and increasing system sizes suggest that 
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FIG. 2: Force distribution for two dimensional systems, (a) Effect of contact number z on P{f) for disordered packings of 
N = 1000 particles. The inset compares the P(/) for a disordered packing with z = 6 and A*' = 1000 and the hexagonal packing 
for A'^ = 2900. (b) The same data as in (a), now plotted as log P{f) vs tends to a straight curve for large z. The inset shows 
that on a smaller range, all curves look Gaussian, (c) Same data as in (a-b), now on a triple log plot. The range in / over 
which P{,f) looks Gaussian grows with contact number z. (d) For fixed small z = 4.5, P{f) appears to approach a Gaussian 
tail for large TV. 




FIG. 3: Force distribution for three dimensional systems, (a) P{,f) for two disordered and a regular fee packing of TV = 500 
particles, (b) Same, now plotted as function of The dashed line corresponds to a hexagonal packing in 2D, which has a 
Gaussian tail - the tail of P{f) for 3D systems is significantly less steep, (c) Same data, now plotted as function of /^'^ - the 
tails for the P{f) of 3D packings are now straight, (d) The change from 2 to 1.7 is also clearly visible in the triple log plot. 
For a range of system sizes and contact numbers, we robustly find that P(/) ~ exp( — (//A)") with an exponent a ~ 1.7 for 3D 
systems — for comparison we also show the Gaussian distribution for the 2D hexagonal packing. Note that for small systems 
and small contact number (A'^ = 250, z = 9.1), finite size deviations, similar to those observed in two dimensions, can be seen. 



the "kink" in the triple log plots becomes less severe for 
large systems (Fig. md)) — our data are not conclusive 
as to whether this kink will disappear for TV — » oo. 

In conclusion, for two dimensional, frictionless systems, 
the ensemble approach yields that have a tail with 
decays at least as fast as a Gaussian. 

Three dimensional packings — We now turn to three 
dimensional systems, which again have been generated 
using molecular dynamics. Similar to what happens in 
two dimensions, Fig. [3^ shows that P{f) decays faster 
than exponentially, and disordered and regular (fee) 
packings have very similar force distributions. However, 
the decay is now slower than Gaussian and much more 
accurately described by P{f) ^ exp(— (//A)") with an 
exponent a = 1.7 ± 0.1, see Fig. [3Kb-d). This exponent 
has been determined from the triple- log plots of Fig. EKd) 
for a range of contact numbers and system sizes, and in 
all cases the slope is close to a = 1.7 over a decade. 

For comparison we have, in Fig. [3Jb-d), also included 
the result for the hexagonal pack, which is seen to de- 
crease significantly more rapidly than the P(/)'s of the 
three dimensional systems. Surprisingly, we thus find 



that the dimensionality of the packing determines the 
nature of the tail of P{f)- 

The effect of shear stress — From experiments on (two 
dimensional) sheared packs of photoelastic grains, it was 
found that the distribution broadened significantly, and 
developed an exponential-like regime in a range up to 
4(/) uM- The ensemble indeed reproduces this quali- 
tative feature for packs under shear. As can be seen 
in Fig. m however, there does not seem to be a simple 
asymptotic decay. This is because the force anisotropy 
induced by the shear stress yields a variation in (/) de- 
pending on the orientation of the contact IT], ll9| . The 
total P{f) becomes a superposition over all orientations, 
of mixed force statistics, and hence lacks a single charac- 
teristic feature. 

Discussion — In none of the cases we investigated, 
P{f) exhibits an exponential tail. The "Boltzmann" 
type arguments based on conservation of total force, 
which was suggested to be responsible for the exponen- 
tial tail [J could in principle have been applicable 
here: indeed, all contact forces in the ensemble are posi- 
tive and add up to a value proportional to the pressure. 
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FIG. 4: Two dimensional disordered system with z — 5.5 
experiencing a shear stress r = axy/crxx (a) While for 

large r, the tail of P{f) viewed over a limited range broadens 
and may appear exponential (inset), the asymptotic decay of 
P{f) for / > 10 in fact becomes steeper (main panel), (b) 
The same point is illustrated in the triple log plots, which also 
show data for r = 0.1 and 0.3. 

This reasoning, however, does not take into account that 
forces have to balance on each grain. Our results under- 
line the importance of these additional constraints, which 
completely alter the properties of P{f)- 

The force distributions obtained in the ensemble are 
consistent with those obtained in most experimental and 
numerical studies. In contrast, experiments that mea- 
sure forces at the boundaries appear to find an exponen- 
tial decay of P(/) 0, H, This remains a crucial issue 
for the understanding of static granular media, since the 
force statistics provides insight to the proper measure to 
weigh the microscopic configurations corresponding to a 
macroscopic experimental protocol At present, one 
still lacks a relation between characteristics of the sys- 
tem (presence of boundaries [2§| , the relative hardness of 
grains and boundaries [29| . friction and nonsphericity, ...) 
and the force distribution. Within the ensemble theory 
it would be possible address the effect of torque balance, 
which is clearly important for real (frictional) systems 
or even for frictionless non spherical grains — both fric- 
tion and nonsphericity can in principle be included in 
the ensemble. Having seen that the normal force balance 
conditions have different effects in two or three dimen- 
sions, it would be very interesting to see whether torque 
balance could yet again change the nature of the tail. 
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